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Broad-tailed force distributions and velocity ordering in a heterogeneous membrane 

model for collective cell migration. 
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Correlated velocity patterns and associated large length-scale transmission of traction forces have 
been observed in collective live cell migration as a response to a "wound" . We argue that a simple 
physical model of a force-driven heterogeneous elastic membrane sliding over a viscous substrate 
can qualitatively explain a few experimentally observed facts; (i) the growth of velocity ordering 
which spreads from the wound boundary to the interior, (ii) the exponential tails of the traction 
force distributions, and (iii) the swirling pattern of velocities in the interior of the tissue. 

PACS numbers: 87.18.Hf 
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I. INTRODUCTION 



The phenomenon of collective ceU migration arises in 
biological processes of morphogenesis, wound heahng, as 
well as cancer growth, and is an active topic of current 
research interest [ll-Q- To understand the basic features 
in collective cell migration as a response to wound heal- 
ing, two-dimensional monolayer patches of Madin-Darby 
canine kidney (MDCK) cells on dcformable substrates 
have been studied in different experiments 0, mWi ■ 

For a physical scientist, there are many interesting as- 
pects that these experiments reveal. The spatially het- 
erogeneous swarming and swirling velocity patterns ex- 
hibited by the cells, studied by particle image velocime- 
try [6|, lZI , are reminiscent of similar pattern formation in 
active nematics and driven granular matter [9|. As time 
passes, a zone of velocity order starting from the wound 
boundary invades the interior of the MDCK tissue [3] , re- 
minding one of phase ordering kinetics [lO| . On the other 
hand, another set of experiments [J, |8| have shown that 
the local traction forces exerted by the MDCK cells on 
the substrate have large fluctuations — the distribution 
of the forces being non-Gaussian with distinct broad ex- 
ponential tails, akin to force distributions in static granu- 
lar piles [il|- Yet, the MDCK cells forming the tissue are 
held to each other and to the substrate by a network of 
Cadhcrin and Intcgrin proteins, respectively |ll. Il2l|. and 
they self-generate active forces due to internal Actin and 
Myosin dynamics. Thus, at the microscopic level, they 
show no resemblance to mechanically driven loose gran- 
ular rods or discs. Needless to say, it is quite a challenge 
to model every observed feature of the MDCK tissue sys- 
tem, as seen in different sets of experiments. In this pa- 
per, we propose a simple statistical-mechanical model for 
the system and show that in can simultaneously give a 
qualitative explanation of the growth of velocity ordering 
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with time, and the large force fluctuations. 



It is well known that the behaviour of large cell col- 
lectives [iJl is qualitatively distinct from that of a single 
cell [i5|. Attempts have been made to model cell as- 
semblies incorporating signal transduction and cell-cell 
signaling in Dictyostelium discoideum [l^|. Collective 
cell migration studies incorporating cell division has been 
done IjJrjSlJ. but the experiments that we are concerned 
with [J, Is lZI have noted, that over the relevant time- 
scales, the growth of cell number via cell division is not 
expected to play a role in the features of interest in this 
paper. The geometrical instabilities such as fingeringand 
tip-splitting of the wound boundary in experiments [a, |7[ 
have been theoretically modelled using ideas of interface 
growth kinetics QiIML The velocity patterns in the inte- 
rior of the cell-sheet [y, 0] have recently been studied by 
a mechanical model, where the cells with a local orienta- 
tion field collectively behave in a viscoelastic fashion [^ . 
We show below that a much simpler model compared to 
these can nevertheless capture two key features of the 
MDCK tissue system, namely, the velocity ordering in- 
vasion, and the large force fluctuations. In particular, 
the novel aspect of our model is to show the crucial role 
that heterogeneity plays in determining the large force 
fluctuations. 



The paper is arranged in the following fashion. In sec- 
tion |TT] we describe our model in details, and justify its 
various components on the basis of relevant experimental 
findings. We also specify the equations and the param- 
eter values used for the simulation study of the model. 
In section IIIIt \-C we describe the three major findings 
of our model. Numerical estimates to make quantita- 
tive connections with the experiments are detailed in sec- 
tion IIIIDI The significance of our model and findings in 
comparison with previous theoretical works are discussed 
in section IIVI 




FIG. 1: A schematic picture (top view) of the deformed dis- 
cretized membrane modeled as a tihed square lattice with 
undeformed spring-length ao of unity. The wound boundary 
cells move towards — ve X. The cells are denoted by the circles 
and the active forces are denoted by arrows (size proportional 
the force magnitude) on them. The varying spring stiffnesses 
are denoted by multiple colors. The overall average substrate 
resistance is denoted by Drag. 



II. MODEL AND PARAMETERS 

We model the cell sheet as an elastic membrane, al- 
though some earlier studies have treated cell monolayers 
as viscoelastic [3, [l^. The motivation for this choice 
comes from the experimental finding [7| that the aver- 
age distance between MDCK cells does not change for 
many hours. It was clearly concluded that the move- 
ments of the cells are correlated showing very limited re- 
arrangement, making the monolayer locally more elastic 
than viscous [7|. Similarly, a recent experimental study 
of internal stresses in cell monolayers [5| have also as- 
sumed the latter to be elastic. In our model, the cells are 
represented as discrete points in continuous space, and 
the cell-cell cadherin connections are represented by sim- 
ple harmonic springs of effective stiffness kq. The elas- 
tic membrane is considered heterogeneous with kq drawn 
from a probability distribution function (p.d.f) P{kq). 
Since there is inherent randomness in the number of in- 
tercellular cadherin connections [221, ^^^ the possibility 
of any such connection to be in multiple distinct struc- 
tural states with different mechanical properties [23| , our 
assumption of heterogeneity seems reasonable. Next, it 
is expected that the connections of the cells with the sub- 
strate via the integrin proteins will break and remake as 
the tissue advances [Ij]. For a single cell, it has been 
theoretically demonstrated that this complicated process 
of cell-substrate interaction can be effectively replaced by 
a linear viscous drag [2al- We extend the latter at the 
tissue level and assume that the substrate exerts a local 
drag force — coV^ where V^ is velocity of the i-th cell 
and Co is the drag constant. The average position of the 
"wound boundary" (see Fig. 1) defines the ^-direction, 
and the direction orthogonal to it will be called X. For 



simulating, we take a. N x N tilted square grid of points 
(Fig. 1), but while thinking of a continuum limit, we will 
assume iV — > oo such that effectively the wound bound- 
ary will be very far from the center of the tissue (as in 
actual experiments 0, 0] ) ■ Finally, we supply the "live 
thrust forces" (originating from the cytoskeletal acto- 
myosin activity in the cells) by hand in two alternate 
ways: (i) Cell— i with position R^ = (Xi,Yi) is given a 
space and time dependent random force F^ = Favo.i + ''7i, 
with Favc.i = ~^o exp(— ni/^)x. Here rii is the row num- 
ber, from the boundary, of the i"^ cell. The boundary row 
is numbered 0, and ^ is a length scale. The components of 
rji = {r]x,r]Y)i are Gaussian white noise with zero mean 
and {r]a.iiti)rij3j{t2)) = 2acx.p{-ni/S,)SijSa.i3S{ti - ^2) 
with a and /3 taking values X, Y. We will refer to 
this as participatory model (PM) (see Fig. 1). Thus 
in the PM model, the magnitude of the noise on force, 
just like the mean force, decays with increasing distance 
from the boundary, (ii) Only the cells at the wound 
boundary row are given a space and time dependent 
random force with an average force Favo,i = —FoSmfi'^ 
added to Gaussian white noise rji, with zero mean and 
{VaA*i)Vfs.j{t2)} = 2aSij5a.pS{ti - t2) (o = {X,Y}, 
/3 = {X, y}). We will refer to this as the leader driven 
model (LDM). 

The motivation for comparing PM versus LDM comes 
from the discussions in Ref. [26[ following the experiment 
of Ref. Q . The mathematical equation used to simulate 
the system is: 
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The inertial term has been dropped as the system is 
clearly overdamped. The index j in the sum in Eq. [T] 
goes over nearest neighbours of i and Kq is the random 
stiffness constant of the spring connecting i and j. The 
unit vector eij ~ (Rj — Ri)/|Rj — R^j, and oq is the 
undeformed length of the spring. 

Here we report numerical results with parameter val- 
ues: Co = 10, N = 128, Fo = 1, cr = 0.2, ao = 1, 
and P{kq) is a uniform box distribution between 0.25 to 
0.75. We will show later that these choices of parameters 
lead to reasonable correspondence to experiments. In the 
presence of the active applied force F^, we time-evolve the 
positions R^ using Eq. [1] and obtain the velocities V^ as 
dRi/dt. The simulation results for Figs. 2, 3, and 4, are 
obtained by providing the cells with zero initial veloci- 
ties, and initial random displacements (from the equilib- 
rium positions) with components uniformly distributed 
over [—5,(5]. This may be expected to be a generic ini- 
tial condition for the cell collective. Periodic boundary 
conditions are assumed along Y. 



III. RESULTS 

In this section we present the three important results 
of our model. 
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FIG. 2: Velocity ordering in the cell-sheet using 5 — 0.025. 
(a) Velocity of the cell points as a function of row number for 
two different times (100 and 1000). The circles and squares 
represent simulation results for PM (^ — 5) and LDM, re- 
spectively, and the fitting lines correspond to the theoretical 
model in Eq. 2. (b) The order parameter (j}„ at spatial row po- 
sition X of the cell at different times t. The a:-range is — 120 
and t-range is — 1000. The velocity ordering front is clearly 
visible, (c) The order parameter $ for the complete tissue, 
as a function of time t, with a fit of t^ " and t for the LDM 
and PM, respectively, (d) Time dependence of experimental 
order parameter in [31 fitted with t and t^". 



A. Invasion of Velocity Ordering 

We first proceed to show how under the action of ac- 
tive forces the cell velocities acquire a bulk ordering. As 
the active forces are preferentially oriented along — x (for 
both the LDM and PM), the velocities of the cells, start- 
ing with the ones close to the boundary and followed 
by the ones in the bulk, gradually orient themselves to- 
wards —X. This velocity ordering is shown in Fig. 2(a) — 
with increasing time, for both the LDM and PM, mean 
row velocity ~Vx of deeper layers (with larger n) rise 
in magnitude. Thus an order-disorder boundary moves 
towards larger n. Interestingly, there is a quantitative 
difference between the two models — the shape of the 
—Vx{n) curve is Gaussian for LDM and has an exponen- 
tial tail for PM. A local order p aram eter for every row 
n can be defined as 0„ = —Vx,i/\'Vi\ (with all velocities 
—Vx,i < 10~^^ set to zero to avoid spurious contribu- 
tions). In Fig. 2(b), 0„ is plotted (with larger mag- 
nitude corresponding to lighter color) as a function of 
space X and time t for LDM — the movement of the 
order-disorder boundary with increasing t is clearly seen. 
A similar plot was observed experimentally in [7|. Next, 
a global order parameter for the whole system can be de- 
fined as $ = {(j>n)n- In Fig. 2(c), $ is shown to increase 
as ~ i^/^ for LDM, and ~ t for PM. To compare with 
the experiments, we have plotted the experimental data 
for $ from [31 in Fig. 2(d); two curve-fits of ^ t^/^ and 
~ t are put against the data, showing that both these 



forms work reasonably well. Thus we have demonstrated 
numerically that our models have similar velocity order- 
ing as seen experimentally in MDCK tissues [3] during 
wound healing. We will now proceed to understand ana- 
lytically which ingredients of our model are essential for 
the above phenomenon, and in particular the reason for 
quantitative differences between LDM and PM. 

It is interesting to note that the growth of Vx (t) above 
can be understood analytically from an analogous 1- 
dimensional problem. We can solve the one-dimensional 
problem of a pulled, non-disordered, elastic chain: 






(2) 



The above equation is a simple 1-D, linearized, contin- 
uum version of Eq. [TJ Here u{x,t) is the displacement of 
any cell, x{€ [0,L]) is the continuum space variable cor- 
responding to the row number ri^, and the non-random 
force is F = —Foexp{—x/^), with du/dx = 0|2;=o (for 
PM), and F ^ Q, with du/dx = Fq/kq\^=q (for LDM). 
The initial condition is taken as u{x, 0) = both for PM 
and LDM. Eq.[2]can be solved for velocity v{x, t) ~ du/dt 
in the limit of large system size {L — >■ oo) with the bound- 
ary condition u{oo,i) = 0. This gives unique analytical 
solutions: 
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scaled dimensionless space and time respectively, and Erf 
and Erfc refer to Error and Complementary Error func- 
tions [23] , respectively. The profile of the ordered veloc- 
ity as a function of x in the 1-D LDM model (Eq. [4]) 
is clearly Gaussian, while in the PM model it is mod- 
ified (Eq. [3]) to have an exponential profile for large 
X. The agreement in mathematical forms between the 
1-D analytical result and the 2-D simulation result in 
Fig. 2a shows that the velocity ordering phenomenon is 
not particularly dependent on stiffness randomness or 
noise, and its essence is captured even in a 1-D prob- 
lem. As can be seen from Fig. 2c, the PM model gives 
a growth law ^ ^ t for a transient period, the reason 
for which may be understood from Eq. 3. For large x, 
ln(— w(a;)) ^ r — x -(- Constant, to leading order, while for 
small X, \n{—v{x)) r^ t — xP' f{T) + Constant (see PM in 
Fig. 2a), where /(t) is a function of r. Thus there is an 
short distance quadratic profile, followed by a long dis- 
tance linear profile in x. Spatial expanse of the quadratic 
profile keeps increasing, and beyond a certain time the 
growth law in PM will become ^ t^/^ just like LDM. 

We would now like to see if a drive from the boundary 
cell layer without any participation from the cells in the 
bulk can invoke other experimentally observed phenom- 
ena. Hence, in the rest of the paper we focus on LDM. 
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FIG. 3: Distributions P(Tx) and P{Ty) of traction forces 
Tx and Ty (in units of 10"^) for LDM, at t = 300 and 
S = 0.025. Gaussian fits are shown for uniform Kq in (a) 
and (b). For random Kq , deviation from Gaussianity, and 
resulting asymptotic exponential tails are shown in (c) and 
(d). The data is for layers 18 (filled symbols) and 19 (empty 
symbols). 



FIG. 4: Distributions P{Tx) and P(Ty) of traction forces 
Tx and Ty (in units of 10"^) for LDM, at t = 300, for 
a distorted lattice with random bond lengths. The uniform 
stiffness constant fco = 1. Gaussian fits are shown for e = 0.15 
in (a) and (b). For e = 0.45, deviation from Gaussianity, and 
resulting asymptotic exponential tails are shown in (c) and 
(d). The data is for layers 18 (filled symbols) and 19 (empty 
symbols). 



B. Traction Force Fluctuations 

A second interesting result of our model is that the 
traction force fluctuations are unusually large as in the 
experiments in Ref. [3]. In our model, the local effec- 
tive traction force on any cell is T^ = F^ — cqV (see 
Eq. [1]). In LDM, F^ = (for any non-boundary cell) 
and so T's are proportional to the local velocities V. 
The probability distributions of the components Tx and 
Ty for the LDM are shown in Fig. 3. For homogeneous 
membrane (constant Kq), the distributions are clearly 
Gaussian (Figs. 3(a-b)). To make sure that the latter 
is not a trivial consequence of the Gaussian distributed 
random thrust forces in the boundary layer, we checked 
the traction distributions when the boundary forces were 
drawn from a (i) box, and (ii) an exponential distribu- 
tion. In both of these distinct cases, we found (data not 
shown) that the traction force components T^ and Ty 
are Gaussian distributed. Thus there is no doubt that 
Central Limit Theorem (CLT) is valid and due to it, the 
local forces in the bulk (being sum of random neighbour- 
ing forces) turn out to be normally distributed. On the 
other hand, for a heterogeneous membrane (random Kq) 
the distributions develop exponential tails (Figs. 3(c-d)), 
indicating a departure from the CLT. The shape of the 
curves of P(Tx) (with mean at T^ ^ 0) and P{Ty) 
(with mean at Ty = 0) have qualitative resemblance to 
experiments — in particular, the widths decrease with 
increase of distance from the wound. 

Breakdown of CLT in [J] is a priori quite intriguing. 

like 
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It was speculated in Ref. [J|, that a g-model 
mechanism maybe at play. Recognizing that the tissue 
is heterogeneous, here we are specifically suggesting that 
an effective g-modcl like mechanism may arise due to un- 
equal stress propagation and accumulation mediated by 



the random cadherin connections. Since the possibility of 
manipulating the strengths of cadherin connections has 
been experimentally demonstrated [5| , our result is open 
to experimental test. 

The membrane can also be made heterogeneous in an- 
other way by introducing variable bond lengths of the 
cell-cell connections, while keeping the spring stiffness 
homogeneous. This naturally leads to deviation from 
regular lattice symmetry considered so far. The bond 
lengths Cg'' are made disordered by imparting new ran- 
dom equilibrium positions to the cells. To do so, the cells 
arc shifted from the regular lattice by e^ (along-X) and 
ey (along- Y), where e^ and ey are drawn from uniform 
distribution over [— e, e]. As can be seen from Eq. 1, the 
magnitude and direction of the spring force is indepen- 
dent of the bond- length (to the first order). Thus for 
"small" lattice distortion the distributions for Tx and 
Ty are expected to be the same as that of the uniform 
non-distorted lattice. This is seen in our simulations for a 
choice of e = 0.15 — the traction distributions are indeed 
Gaussian (see Figs. 4(a) and 4(b)). Contrary to this, with 
"large" random lattice distortion, one would expect from 
Eq. 1 random harmonic forces, leading to a departure 
from the latter result. We indeed see this when we make 
e large (say 0.45) — the traction distributions develop ex- 
ponential tails as shown in Figs. 4(c) and 4(d). Thus in 
two types of membrane heterogeneity — random spring 
stiffnesses and random bond lengths — we have found 
that non-Gaussian traction force distributions arise. 



C. Swirling Patterns in the Bulk 

We now turn our attention to velocity patterns which 
develop in the bulk of the system. The cells in the con- 
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FIG. 5: The velocity-velocity correlation function Cw against 
scaled distance n along y, at f = 500 (•) and t = 1000 (o) 
for LDM using S = 0.45. Inset: Velocities of cells in the bulk 
region (dimension 8x8) are shown (at t = 1000) — swirls 
can be clearly seen. 



fined region of the experimental setup in rcfs. [6|, lZ| are 
expected to be under internal stress due to the confine- 
ment from the boundary before the cell sheet is allowed to 
expand. In order to mimic this internal stress, we provide 
generic random initial positions to our cell lattice. The 
elastic relaxation of this "pre-strained" sheet help ex- 
cite spatial modes through the harmonic couplings. The 
smaller wavelength modes damp out faster, leaving large 
wavelength velocity swirls at late times (inset of Fig. 5 
for LDM). To precisely quantify the correlations in these 
patterns, we show the velocity-velocity correlation func- 
tion Cvy{n) = (v(no).v(no + n))/(v^) in Fig. 5; here 
n = y/-\/2ao is the scaled distance in units of average 
inter-cellular spacing along Y. The average (• • • ) is done 
over ensembles as well as cell locations no, belonging to 
a strip in the bulk region where (v(no)) = 0. We note 
that Cvv (n) shows change of sign beyond some layers (as 
seen also in experiments [231) reflecting the bending of 
velocity field over space. The correlation range increases 
with the time of observation as expected — it is ss 8 — 10 
cell layers which is roughly similar as distances seen in 
experiments [3, [2^. The correlations observed in ref. 
[28| are certainly influenced by substrate deformability 
and cell birth. But the fact that such correlations are 
also observed otherwise [3, 01 indicates that they may 
not be the only factors governing the swirling patterns. 
Although our model cannot make these distinctions we 
note that it elucidates the role of inherent elasticity and 
pre-strain of the cell-sheet in producing such patterns. 



D. Numerical Estimates of Length, Time and 
Force Scales 



We have shown so far that a boundary layer driven het- 
erogeneous elastic sheet can produce qualitatively many 
experimental observations in collectively migrating ep- 
ithelial cells. We now make numerical estimates of vari- 
ous quantities to find out whether our results are quanti- 
tatively meaningful in comparison with the experiments. 



As shown in Fig. 1, the bond- length ap is unity. In real 
units it may be taken as 20yLtm (the average cell-cell sepa- 
ration Q). To estimate the unit of time to, we start from 
Figs. 2(c) and (d). For the LDM in 2(c), the slope of 

— 1/2 

the line in units of t^ is 0.0139, while the experimen- 
tal data in 2(d) has a slope 0.139 hr~^". Equating the 
two, gives to w O.Olhr = 36s. This tells us that the non- 
dimensional velocity in Fig. 2(a) in the boundary layer 
(for LDM) - 10~2 corresponds to IQ-'^ao/to « 20Aim/hr. 
This is in the ball-park of the velocities quoted in exper- 
iments Q. The curves in Fig. 2(a) are for times t = 1 
hr (100 in simulation) and 10 hrs (1000 in simulation). 
To make contact of forces in Fig. 3 with experiments, we 
choose to first obtain the bond-stiffness kq in real units. 
The properties of a continuum cell sheet may be chosen 
as those appearing in the supplementary material of [5|: 
Young's modulus E ^ 10 kPa, Poisson's ratio ly = 0.5, 
and sheet thickness of /i = 5/im. We first consider a cell 
sheet of dimension oq x ao x /i, on which homogeneous 
tension 1 Pa is applied. The relative change in surface 
area AA/A = 2{l-u)/E = 10"". On the other hand, ap- 
plying equivalent edge force of F = IPa x oq x ft, = 10^^° 
N on a flo X ao square, whose sides are made of springs of 
stiffness kq, we get AA/A = F/aoKo- Equating the two 
relative area changes, we get kq = 0.05 N/m. But we 
have used an average stiffness value of 0.5 in dimension- 
less units in our simulation. This gives the actual force 
unit to be /o = 0.05 x ao/0.5 N = 2 x 10"*^ N. In dimcn- 
sionless units our forces in Fig. 3 for the 19"^ layer are 
~ 0.5 X 10~^, which is equivalent to 10~^ N. This implies 
a traction of 10~^/ao = 2.5 Pa, which is one order of 
magnitude lower than that reported in [^. The dimen- 
sionless time 300 reported in Fig. 3 translates to sa 3 
hrs. Using ao, to and /o, we see that the dimensionless 
value of Co = 10 used in our simulation, is equivalent to 
10 X foto/ao ~ 36 N-s/m in real units. This can be con- 
verted into drag co-efficient C — co/o-o — 25 pN-hr//im'^, 
which is an order of magnitude lower than reported in 
Q. In Fig. 5, at our simulation time t — 1000 (equiva- 
lently « 10 hrs) the length scale associated with the first 
minimum of the correlation curve is 5 x ao = lOO/^tm. In 
[231 at around 10 hrs, a similar reported length scale is 
200//m (or « 10 cell layers), which is higher than our re- 
sult only by a factor of two. Thus we see that, although 
our model is simple, it can make reasonably close contact 
to experiment even quantitatively. 

One serious drawback is that our boundary tractions 
are very large ^ lO'^ times compared to the bulk — in 
reality Q forces do not diminish so fast. It would be 
interesting to modify LDM in the future by incorporating 
cellular thrusts from the bulk, to see if the description 
becomes more realistic in a quantitative sense. 



IV. DISCUSSION AND CONCLUSION 

Recent experiments on collective migration of MDCK 
cells have thrown up several interesting puzzles, which in 



turn have spurred various theoretical modeUng attempts. 
Before we summarize the main results of this paper, we 
would like to situate our work with respect to the con- 
tributions made by the earlier theoretical models. The 
model of Ref. [3| treats the cell sheet as a viscoelastic 
medium supplemented with a director field to describe 
the local cellular orientations. This model obtains the 
dependence of the velocity of the wound boundary on the 
viscoelastic parameters of the cell-sheet, and also shows 
complex correlated velocity patterns in the bulk. An- 
other model |2| concentrates specifically on the dynamics 
of the boundary of the cell-sheet. By introducing a com- 
petition between the curvature dependent driving force, 
and the elastic and viscous resistance of the cell sheet, 
the fingering instability as seen in experiment [7| is re- 
produced by this model. In contrast to these models, we 
treat the cell-sheet as an elastic membrane, and hope to 
capture some of the phenomena at early times. This is 
motivated by a direct experimental observation [7|, and 
supported by treatment of cell-sheet as an elastic mate- 
rial in another set of experiments [5| . We note that the 
phenomena we address in this paper, namely, growth of 
bulk velocity order parameter [7[ , and traction force dis- 
tributions y have not been addressed in the aforemen- 
tioned publications 0, Hi ■ At the same time, the LDM 
cannot produce pronounced fingering due to the lack of 
flowy behaviour in our model. A very recent paper [20[, 
which studies the effect of cell proliferation and migra- 
tion leading to contact inhibition, introduces a simple 
one-dimensional model, where the cells plastically spread 
in presence of cellular thrust forces from the boundary 
(similar to LDM in our paper). Nevertheless, since this 
model is one-dimensional, quite naturally, it cannot cap- 
ture the two-dimensional phenomenology. 

In this paper we have identified few minimal mechani- 
cal ingredients — heterogeneous elastic membrane, fluid- 



viscous drag, and the active drive of cells from the bound- 
ary to mechanically pull the system — which can explain 
three aspects of collective cell migration: (a) macroscopic 
velocity ordering, (b) breakdown of CLT for traction 
force fiuctuations, and (c) velocity correlations associ- 
ated with swirls. Perhaps the most interesting result is, 
that membrane heterogeneity has the capacity to induce 
broad tails in the distribution of traction forces. The 
mechanism that we propose here is very reminiscent of 
the q- model for static granular assemblies [ll|. At the 
same time, we would like to point out that there are sig- 
nificant differences of our model from the q-model. We 
have a mobile network of cells, velocity dependent dissi- 
pative forces, and tensile force transmissions, as opposed 
to the static transmission of compressive forces in the 
granular assemblies. These differences may invite fur- 
ther analytical exploration of our current model in the 
future. 

The three results in the paper show a close qualita- 
tive resemblance to the experiments. Even the quan- 
titative estimates seem reasonable, albeit with a major 
drawback that the values of traction forces and veloci- 
ties diminish much faster than the experimental values. 
This can be attributed to the fact that our model does 
not pump energy in the interior of the cell layer through 
active cellular thrusts. Recent experiments |5| hint that 
cellular polarizations and cellular active forces are possi- 
bly tied to mechanical stress cues from surrounding cells. 
This demands our model to go beyond being purely me- 
chanical, by including a coupled dynamics (a cross-talk) 
between active cellular forces and mechanical harmonic 
forces. While we would explore these in future, we con- 
clude by noting that this work sets a benchmark by show- 
ing the achievements and limitations of a rather simple 
mechanical model for collective cell migration. 
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